Phase transitions in systems of self-propelled agents and related network models 
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An important characteristic of flocks of birds, school of fish, and many similar assemblies of self- 
propelled particles is the emergence of states of collective order in which the particles move in the 
same direction. When noise is added into the system, the onset of such collective order occurs 
through a dynamical phase transition controlled by the noise intensity. While originally thought to 
be continuous, the phase transition has been claimed to be discontinuous on the basis of recently 
reported numerical evidence. We address this issue by analyzing two representative network models 
closely related to systems of self-propelled particles. We present analytical as well as numerical 
results showing that the nature of the phase transition depends crucially on the way in which noise 
is introduced into the system. 
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The collective motion of a group of autonomous par- 
ticles is a subject of intense research that has potential 
applications in biology, physics and engineering [H, 0, Q . 
One of the most remarkable characteristics of systems 
such as a flock of birds, a school of fish or a swarm of 
locusts, is the emergence of ordered states in which the 
particles move in the same direction, in spite of the fact 
that the interactions between the particles are (presum- 
ably) of short range. Given that these systems are gener- 
ally out of equilibrium, the emergence of ordered states 
cannot be accounted for by the standard theorems in sta- 
tistical mechanics that explain the existence of ordered 
states in equilibrium systems typified by ferromagnets. 

A particularly simple model to describe the collective 
motion of a group of self-propelled particles was proposed 
by Vicsek et al. [J| . In this model each particle tends to 
move in the average direction of motion of its neighbors 
while being simultaneously subjected to noise. As the 
amplitude of the noise increases the system undergoes a 
phase transition from an ordered state in which the par- 
ticles move collectively in the same direction, to a disor- 
dered state in which the particles move independently in 
random directions. This phase transition was originally 
thought to be of second order. However, due to a lack of 
a general formalism to analyze the collective dynamics of 
the Vicsek model, the nature of the phase transition (i.e. 
whether it is second or first order) has been brought into 
question 

In this letter we show that the nature of the phase 
transition can depend strongly on the way in which the 
noise is introduced into these systems. We illustrate this 
by presenting analytical results on two different network 
systems that are closely related to the self-propelled par- 
ticle models. We show that in these two network models 
the phase transition switches from second to first order 
when the way in which the noise is introduced changes 



from the one presented in [1] to the one described in 

The first network model, which we will refer to as 
the vectorial network model, consists of a network of 
N 2D-vectors (represented as complex numbers), {<j\ = 
e Wl ,a 2 = e i02 ,...,a N = e i0N }, all of the same length 
|er n | = v and whose angles {0i(t), O^it), . . . , 0jv(£)} can 
change in time. Each vector o~ n interacts with a fixed 
set of K other vectors, {cr ni , . . . , a nK }, randomly chosen 
from anywhere in the system. We will call this set of K 
vectors the inputs of o~ n . Once each vector a n has been 
provided with a fixed set of K input connections, the dy- 
namics of the network are then given by one of the two 
following interaction rules: 

9 n (t + 1) = Angle J -L £ (f) 1 + ^( t ), (1) 

9 n (t + 1) = Angle I -L J2 (*) + »?e* W J , (2) 

where for any vector v = \v\e lc ^ we define the function 
Angle (i/) = 4>, and £(t) is a random variable uniformly 
distributed in the interval [— tt, 7t]. The the dynamics of 
the network is fully deterministic for rj = and becomes 
more random as the parameter r\ increases. In what fol- 
lows, we will refer to the quantity (1 / vK) Y)j—y Cn^ as 
the average contribution of the inputs of o~ n . 

To quantify the amount of order in the system we de- 
fine the instantaneous order parameter ip(t) as 



ib(t) = lim — — 
JV^oo vN 



X>n(t) 



(3) 



In the limit t — ► 00, the instantaneous order parameter 
tp(t) reaches a stationary value ip [3, 0, IS 0] ■ Thus, in 



the stationary state all the vectors are aligned if if) ~ 1, 
whereas if tp ~ the vectors point in random directions. 

The interaction rules given in Eqs. (p} and ([2]) were 
proposed by Vicsek et al. in Ref. [4[, and by Gregoire 
and Chate in Ref. [H, respectively. The difference be- 
tween these two interaction rules consists in the way in 
which the noise is introduced: In Eq. (fTJ) the noise is 
added outside the Angle function, i.e. after the Angle 
function has been applied to the average contribution of 
the inputs. On the other hand, in Eq. ([2]) the noise is 
added inside the Angle function, i.e. it is added directly 
to the average contribution of the inputs. In Ref. [5j, 
Gregoire and Chate posed the question as to whether 
these two rules lead to the same type of phase transition. 

In this letter we show that the interaction rules in 
Eq. {1} and Eq. ([2|) produce different types of phase tran- 
sitions in the network systems under consideration, which 
suggests that a similar effect is being observed in [j| for 
the self-propelled systems. 

Obviously, in the self-propelled particle models the el- 
ements do not interact through a network. Instead, they 
move in a 2D space, each particle interacting locally with 
the particles that fall within a certain radius. This mo- 
tion allows particles that are initially far apart to meet, 
interact, and separate again, giving rise to effective long 
range interactions. On the other hand, in our vectorial 
network model the particles are fixed to the nodes of a 
network. The long-range correlations produced by the 
motion of the particles in the self-propelled models are 
proxied in our network model through randomly choosing 
the inputs of each element from anywhere in the network. 
An underlying assumption of our work is that the exis- 
tence and nature of the phase transition depends mostly 
on the occurrence of such long-range interactions, and 
less crucially on whether they are produced by the mo- 
tion of the particles or by the network topology 0, 0]- 
While the exact relation between these two ways of es- 
tablishing long-range interactions is not yet known, it has 
been shown that a strong parallel can be established be- 
tween them 0, 0. Further, below we show that there 
are at least two limits in which they are fully equivalent: 
for large particle speeds and for high densities (see Fig. [I] 
and Fig. EL 

In Ref. [7| it was proven that, as the noise amplitude rj 
increases, the vectorial network model with the interac- 
tion rule given as in Eq. {1} undergoes a continuous phase 
transition from ordered states where ip > 0, to disordered 
states where -0 = 0. Fig. Q] shows this phase transition 
obtained numerically for N = 20000 and K = 5. It also 
displays the phase transition in the Vicsek model for a 
system with the same N, a density such that the average 
number of interactions per particle is also K — 5, and 
increasing particle speeds. As can be seen from Fig. [TJ 
the Vicsek model curves approach continuously the net- 
work model curve as » — > oo. This supports the idea 
that in both cases a second order phase transition is ob- 
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FIG. 1: (Color online.) Phase diagram of the Vicsek model 
and the vectorial network model for the case in which the noise 
is added as in Eq. JTJ . When the speed v of the particles in the 
Vicsek model increases, the phase transition converges to that 
of the vectorial network model. The numerical simulations 
were carried out for systems with N = 20000 particles and an 
average number of interactions per particle K = 5. 



served when the noise is introduced as in Eq. |T]), albeit 
the finite size effects observed near the critical point. 

The probability distribution function (PDF) of the 
sum ^ Si=i a rij (t) + rye^O that appears in Eq. ([2]) is 
computed as for a random walk assuming that all the 
terms are statistically independent. By projecting this 
PDF onto the unit circle we can establish a recursion re- 
lation for the order parameter, which for K ^> 1 becomes 
+ 1) = M v (ip(t)), where 



Mn(m) = 



^^(1,1,2,!^) if^)<, 

2^1 (Vii,^) ifr,<m 



(4) 

and 2-F1 (a, b; c, x) is the Gauss hypergeometric function. 
The function M v (ip) is shown in Fig. [2^ for different val- 
ues of r\ (solid curves). This figure also displays with 
symbols the numerical dynamical mapping computed for 
the self-propelled model with the interaction rule given 
in Eq. ([2]), V = 20000 particles, and an average num- 
ber of interactions per particle K — 100. Clearly, the 
numerical mapping coincides with the theoretical result 
for Ai v (ip), showing that the network and self-propelled 
systems are also equivalent in the high density limit case 
considered here. The numerical mappings for the self- 
propelled system were obtained by placing the particles 
in various random initial conditions constrained to pro- 
duce every order parameter value ip(t) in the x-axis, and 
then computing one time step using Eq. ([2]) to obtain the 
corresponding value ip(t + 1) in the y-axis. 

The fixed points of the dynamical mapping tp(t + 1) = 
M.ri(tp(t)) give the stationary values of the order param- 
eter. From Eq. ([4]) it is clear that ip — is always a 
fixed point. However, the stability of this fixed point 
changes depending upon the value of r\. By numerically 
solving Eq.(|3|) to obtain the fixed point, we find that 
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FIG. 2: (Color online.) Graph of the dynamical mapping 
M v (tp)- (a) The interaction rule is as in Eq. ([2]). The solid 
curves correspond to the analytical solution given in Eq. Q, 
and the symbols to the numerical simulation carried out for a 
system with N = 20000 and an average number of interactions 
per particle K — 100. (b) The interaction rule is as in Eq. fTJ. 
The curves were computed numerically for a system with N = 
20000 particles and K = 5. In (a) the non-zero stable fixed 
point appears discontinuously as rj decreases, whereas in (b) 
it appears continuously. 



for 0.672 < rj the only stable fixed point is ip = 0. As 
7] — > 0.672 from above, the graph of M. v (ip) moves closer 
to the identity and eventually another non zero stable 
fixed point ip' appears discontinuously when rj « 0.672 
(see the point indicated with an arrow in Fig. [2^). For 
1/2 < rj < 0.672 there are actually two stable fixed 
points. In this region of bi-stability the system shows 
hysteresis. Finally, when rj < 1/2 the fixed point ip = 
becomes unstable and only the non zero fixed point ip' 
remains stable. Contrary to this, when the noise is in- 
troduced as in Eq. ([1]) the non-zero stable fixed point 
appears continuously, as can be seen from Fig. [2b. 

The validity of these results is corroborated by numer- 
ical simulations carried out for networks with N = 10 5 
and K = 20. Fig.[3]shows the fixed point ip of Eq. ((4]) as a 
function of the noise intensity rj (solid line) . The disconti- 
nuity of the order parameter ip at rj = 0.672 and rj = 1/2 
is apparent. The dashed and dotted-dashed curves are 
the plots of the results from the numerical simulation for 
the cases in which all the vectors were initially aligned in 
the same direction (ip(0) = 1), and when the vectors were 
initially oriented in random directions (ip(0) — 0), respec- 
tively. In the region of bi-stability 1/2 < rj < 0.672, the 
system reaches one or the other of the two stable fixed 
points depending upon the initial condition. 

The theoretical curves presented in Fig. [3] show the 
"limits of metastability" for the system; i.e. the maxi- 
mal and minimal values of rj for which the system has 
bistable behavior (hysteresis). Clearly, specific realiza- 
tions of the system cannot be driven all the way to the 




FIG. 3: (Color online.) Phase diagram of the vectorial net- 
work model for the case in which the noise is added as in 
Eq. (T2J). The solid line corresponds to the prediction obtained 
from Eq. Q. The dashed and dotted-dashed curves are the 
results of the numerical simulation starting out the dynam- 
ics from initial conditions for which ip(0) w 1 and ip(0) « 0, 
respectively. The phase transition in this case is clearly dis- 
continuous. 



limits of metastability and decay at values of 77 slightly 
above 1/2 and below 0.672, as observed in the graph. 
Additionally, Eq. (01 was obtained in the limit of large 
K, however, already for K = 20 the agreement is good. 

The second model that we consider is a majority voter 
model in which the network elements o~ n can acquire only 
two values, +1 or —1. We can think of this system as 
a society in which every individual o~ n has to make a 
decision about an issue with two possible alternatives, 
either +1 or — 1. Again, each element a n receives inputs 
from a set of K other elements randomly chosen from 
anywhere in the system. Let us first consider the case 
in which the interaction between a n and its K inputs, 
}, is given by 



<7„(i + l) = Sign 



(5) 



where Sign[x] = 1 if x > 0, Sign[x] = — 1 if x < 0, rj is 
a parameter that takes a constant value in the interval 
[0, 1], and £(t) is a random variable uniformly distributed 
between [—1,1]. For the Sign function to be well defined 
we choose K as an odd integer. Eq. ([5]) is similar to 
Eq. |T]) in that the noise is added to the sign of the aver- 
age contribution of the inputs. Since in this case o~ n is a 
discrete variable that takes only the two values +1 or -1, 
the Sign function has to be applied again. This interac- 
tion rule reflects the fact that an individual in a society 
usually tends to be of the same opinion as the majority 
of his "friends" (inputs), though, with probability rj/2 he 
can have the opposite opinion. 

The instantaneous order parameter ip(t) is defined as in 
Eqs. but now the vertical bars represent the absolute 
value instead of the norm of a vector. 

In Ref. [H it has been shown that the majority voter 
model with the interaction rule given by Eq. ([5]) under- 
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FIG. 4: (Color online.) Phase transitions in the majority 
voter model when the noise is added as in Eq. ([5]) (solid line 
with circles), and as in Eq.© (dashed line with diamonds). 
The numerical simulations were carried out for systems with 
N = 10 5 and K = 3. 



goes a continuous phase transition as the value of rj is 
varied. In Fig. [4] we reproduce this phase transition for 
networks with N = 10 5 and K = 3 (solid curve with cir- 
cles). It is clear from this figure that, when the noise is 
added as in Eq. I|5|). the phase transition in the majority 
voter model is indeed continuous. 

Let us now consider the case in which the interaction 
rule between a n and its inputs is given by 



a„(t+l) = Sign 



1 K 



(0) 



where £(t) is a random variable uniformly distributed in 
the interval [—1,1] and r\ G [0,1]. This rule is similar 
to that given in Eq. @ in that the noise is added to the 
average contribution of the inputs and then the Sign func- 
tion is evaluated. Again, the PDF of the sum appearing 
in Eq. ^ can be computed as for a random walk as- 
suming that all the terms are statistically independent. 
Integrating the PDF over all positive values of the sum 
we obtain a recursion relation for the order parameter 
(see 0), which for K = 3 becomes 



for < r\ < | 
fOT \ < V < \ 



for | < r\ 



Although ip = is always a fixed point of the previous 
equation, a stability analysis reveals that for < r] < 1/2 
the solution ip = is unstable and the only stable fixed 
point is ip = 1. In the region 1/2 < T) the fixed point 
= 1 disappears altogether and the only fixed point 
that remains is ip = 0, and it is stable. Therefore, at 
the critical value 77 = 1/2 the system undergoes a dis- 
continuous phase transition from a totally ordered state 
(ip = 1) to a fully random state (ip = 0). Fig. (j4|) shows 
the results of the numerical simulation for a system with 
N = 10 5 and K — 3 (dashed curve with diamonds). 



In summary, we have analyzed numerically and an- 
alytically the phase transition from ordered to disor- 
dered states in two network models that capture some 
of the main aspects of the interactions in systems of 
self-propelled particles. In particular, the self-propelled 
model becomes equivalent to the vectorial network model 
in the limit of large speeds or high densities. We have 
shown that for the two network models, the phase tran- 
sition changes from second order to first order depending 
on the way in which the noise is introduced into the sys- 
tem. This change is consistent with the results reported 
by Vicsek et al. in [4], and by Gregoire and Chate in 
[5(. This consistency suggests that a similar effect is be- 
ing observed in the self-propelled model and motivates 
a deeper analysis in order to determine the nature of its 
phase transition. Clearly, the two ways of introducing 
noise correspond to different physical situations. On the 
one hand, with the Vicsek type of noise the uncertainty 
falls on the decision mechanism. On the other, introduc- 
ing the noise a la Gregoire-Chate, the decision function 
is perfectly determined and the uncertainty falls on the 
the arguments of this function. There is no reason to ex- 
pect, a priori, similar behaviors under these two different 
physical situations. 

Note added in proof: While this manuscript was being 
reviewed, Ref. |l0j appeared, in which it is shown that 
the first-order phase transition found in Ref. || by means 
of the Binder cumulant is a numerical artifact. 
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